#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double **aloca_matriz(double **a, int n)
{
    int i;
    a = (double**) malloc(n*sizeof(double*));
    for(i = 0; i < n; i++)
    {
        a[i] = (double*) malloc(n*sizeof(double));
    }
    return a;
}

double *aloca_vetor(double *vet, int n)
{
    vet = (double*) malloc(n*sizeof(double));
    return vet;
}

void imprimi_matriz(double **mat, int n, FILE *f)
{
    int i,j;
    for(i = 0; i < n; i++)
    {
        for(j = 0; j < n; j++)
        {
            fprintf(f,"%lf \t", mat[i][j]);
        }
        fprintf(f,"\n");
    }
}

void imprimi_vetor(double *vet, int n, FILE *f)
{
    int i,j;
    fprintf(f,"\n\n");
    for(i = 0; i < n; i++)
    {
        fprintf(f,"%lf\n", vet[i]);
    }
}
void LU(double **a, int n, FILE *f)
{
    int i,j,k;
    double valor;
    double **U;
    double **L;
    U = aloca_matriz(U,n);
    L = aloca_matriz(L,n);
    U = a;
    for(k = 0; k < n; k++)
    {
        for(j = k; j < n; j++)
        {
            L[j][k] = U[j][k] / U[k][k];
        }

        for(i = k+1; i < n; i++)
        {
            valor = U[i][k] / U[k][k];
            for(j = k; j < n; j++)
            {
                U[i][j] = -(valor * U[k][j]) + U[i][j];
            }
        }
    }
    fprintf(f,"\nL:\n");
    imprimi_matriz(L,n,f);
    fprintf(f,"\nU:\n");
    imprimi_matriz(U,n,f);

}

double *fatoracaoLU(double **a, double *b, int n)
{
    double **L;
    double **U;
    double *x;
    double *y;

    x = aloca_vetor(x,n);
    y = aloca_vetor(y,n);
    int i,j,k;
    double soma;
    L = aloca_matriz(L,n);
    U = aloca_matriz(U,n);

    for(i = 0; i < n; i++)
    {
        for(j = 0; j < n; j++)
        {
            L[i][j] = 0;
            U[i][j] = 0;
        }
    }

    for(j = 0; j < n; j++)
        U[0][j] = a[0][j];

    for(i = 0; i < n; i++)
        L[i][0] = a[i][0]/U[0][0];

    for(i = 0; i < n; i++)
    {
        for(j = 0; j < i+1; j++)
        {
            soma = 0.0;
            for(k = 0; k < j; k++)
            {
                soma =  soma + (L[i][k] * U[k][j]);
            }
            L[i][j] = a[i][j] - soma;
        }

        for(j = i; j < n; j++)
        {
            soma = 0.0;
            for(k = 0; k < i; k++)
            {
                soma = soma + (L[i][k] * U[k][j]);
            }
            U[i][j] = (a[i][j] - soma)/L[i][i];
        }
    }

    for(i = 0; i < n; i++)
    {
        x[i] = 0;
        y[i] = 0;
    }
    y[0] = b[0]/L[0][0];

    for(i = 1; i < n; i++)
    {
        soma = 0;
        for(j = 0; j < i; j++)
        {
            soma = soma + (L[i][j] * y[j]);
        }
        y[i] = (b[i] - soma)/L[i][i];
    }
    x[n-1] = y[n-1]/U[n-1][n-1];

    for(i = n-1; i > -1; i--)
    {
        soma = y[i];
        for(j = i+1; j < n; j++)
        {
            soma = soma - (U[i][j] * x[j]);
        }
        x[i] = soma/U[i][i];
    }

    //printf("\nx:\n");
    //for(i = 0; i < n; i++)
    //{
    //    printf("%lf\n",x[i]);
    //}
    return x;
}

void ler_matriz(double **a, int n)
{
    int i,j;

    for(i = 0; i < n; i++)
    {
        for(j = 0; j < n; j++)
        {
            if(i==j)
                a[i][j] = -2.0;
            else if(i==j+1 || j==i+1)
                a[i][j] = 1;
            else
                a[i][j] = 0;
        }
    }
}

void ler_vetor(double *vet, int n, double h)
{
    int i;
    double x = 0;

    for(i = 0; i < n; i++)
    {
        x = x+h;
        vet[i] = pow(h,2)*(-50*(exp(-30*(pow(0.5 - x,2)))));
    }
}

int main()
{
    FILE *f = fopen("saida.txt","w+");
    double *b;
    double **a;
    double *x;
    double h = 0.125;
    int n = 1/h;

    a = aloca_matriz(a,n);
    b = aloca_vetor(b,n);
    x = aloca_vetor(x,n);

    ler_matriz(a,n);
    ler_vetor(b,n,h);
    x = fatoracaoLU(a,b,n);
    //x = fatoracaoLU(a,x,n);
    LU(a,n,f);

    imprimi_vetor(x,n,f);
    return 0;
}
